include("geometry.jl")

function RHS(Q, Op, Geo, t)
    """
    Purpose: Evaluate RHS in 2D Euler equations, discretized on weak form
    with a numerical flux
    """
    # set basic informations
    K = Geo.K                   # number of elements
    Np = Op.Np                  # number of soulution points at each element
    Nfp = Geo.Nfp               # number of flux points at each face
    Nfaces = Geo.NFaces         # number of faces at each element

    # solution points in interior and exterior elements
    vmapM = reshape(Geo.vmapM, Nfp*Nfaces, K)
    vmapP = reshape(Geo.vmapP, Nfp*Nfaces, K)

    # compute volume contributions
    u = Q
    a = 1; b = 1
    F = a*u; G = b*u

    # compute weak derivatives
    rhsQ = zeros(Float64, size(Q))

    # calculate divergence of the discontinuous flux part

    @inbounds dFdr = Op.Drw*F
    @inbounds dFds = Op.Dsw*F
    @inbounds dGdr = Op.Drw*G 
    @inbounds dGds = Op.Dsw*G
    @inbounds rhsQ = (Geo.rx .* dFdr + Geo.sx .* dFds) + (Geo.ry .* dGdr + Geo.sy .* dGds)


    # evaluate '-' and '+' traces of conservative variables
    QM = zeros(Float64, Nfp*Nfaces, K)
    QP = zeros(Float64, Nfp*Nfaces, K)
    @inbounds QM = Q[vmapM]   # get interior value
    @inbounds QP = Q[vmapP]   # get exterior value

    # set boundary conditions by modifying positive traces
    x = Geo.SPs[:,:,1]
    y = Geo.SPs[:,:,2]    
    QB = exp.(-1*((x .- a*t .- 5.0).^2 + (y .- b*t).^2))
    QP = QB[vmapP]

    # evaluate primitive variables & flux functions at '-' and '+' traces
    fM = a*QM
    gM = b*QM
    fP = a*QP; gP = b*QP

    # compute the Rusonov numerical fluxes
    λ = sqrt(a^2 + b^2)

    # Lifting fluxes
    nx = Geo.Normal[:,:,1]
    ny = Geo.Normal[:,:,2]
    
    numerical_flux = nx.*(fP + fM) + ny.*(gP + gM) + λ*(QM - QP)
    
    rhsQ = rhsQ - Geo.LIFT*(Geo.Fscale.*numerical_flux/2)
    rhsQ
end